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Abstract. A Toeplitz matrix is one in which the matrix elements are constant along 
diagonals. The Fisher-Hartwig matrices are much-studied singular matrices in the 
Toeplitz family. The matrices are defined for all orders, N. They are parametrized by 
two constants, a and (3. Their spectrum of eigenvalues has a simple asymptotic form 
in the limit as N goes to infinity. Here we study the structure of their eigenvalues 
and eigenvectors in this limiting case. We specialize to the case with real a and and 
< a < \(3\ < 1, where the behavior is particularly simple. 

The eigenvalues are labeled by an index I which varies from to N — 1. An 
asymptotic analysis using Wiener-Hopf methods indicates that for large N, the jth 
component of the Ith eigenvector varies roughly in the fashion lni/)j ip l j + 0(1/N). 
The Ith wave vector, p l , varies as 

(I) p l = 2nl/N + i(2a + l)\nN/N + 0(l/N) 

for negative values of j3 and values of l/(N — 1) not too close to zero or one. 
Correspondingly the Ith eigenvalue is given by 

(II) e l = a(exp(-V)) + o(l/ N) 

where a is the Fourier transform (also called the symbol) of the Toeplitz matrix. 

Note that p l has a small positive imaginary part. For values of j/N not too close 
to zero or one, this imaginary part acts to produce an eigenfunction which decays 
exponentially as j/N increases. Thus, the eigenfunction appears similar to a bound 
state, attached to a wall at j — 0. Near j = this decay is modified by a set of bumps, 
probably not universal in character. For j/N above 0.6 the eigenfunction begins to 
oscillate in magnitude and shows deviations from the exponential behavior. 

The case of0<a</3<l need not be studied separately. It can be obtained from 
the previous one by a "conjugacy" transformation which takes ipj into i^N-j-i- This 
"conjugacy" produces interesting orthonormality relations for the eigenfunctions. 
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1. Introduction 

1.1. Definition 

A N x N Toeplitz matrix is of the form 





T(0) 


T(-l) 


T(-2) 


T(-3) 




T(l) 


T(0) 


T(-l) 


T(-2) 


T = 


T(2) 


T(l) 


T(0) 


T(-l) 




T(3) 


T(2) 


T(l) 


T(0) 



L ' ' ' J NxN 

A family of these matrices can be generated by doing a Fourier transformation of a 
function a(z) defined on the unit circle 

The generating function a(z) is known as the symbol. 

A Toeplitz operator is a Toeplitz matrix in which N, the number of rows and 
columns is taken to be infinite. As N goes to infinity, the eigenvalues and eigenfunctions 
of the Toeplitz matrix might be expected to converge to those of the corresponding 
operator. However, the convergence is quite non-uniform and subtle as one can see from 
the literature [H [2] . 



1.2. History 

Toeplitz matrices have many applications in physics [3]. One example is that in the 
two-dimensional Ising model, the spin-spin correlation function of the square lattice can 
be written as a Toeplitz determinant [H \5\ E] . More generally they arise whenever a 
line of "impurities" exists in an otherwise uniform system. This line is then represented 
by a matrix in which the interaction between different impurities depends only on the 
distance between them. 

There has been a considerable study of the behavior of these matrices in the limit 
as their order goes to infinity. Szego [7] found an asymptotic determination of Toeplitz 
determinants, which was then extended to a wider class of models by Hartwig and Fisher 
[HE]. The Szego class is defined by symbols which are non-singular on the unit circle. 
These probably have no singularities in their eigenfunctions. The more interesting class 
defined by Hartwig and Fisher has a symbol of the form: 

a Q ^z) = (2-z-l/z) a (-zf (3) 

They have singularities involving irrational powers of N in their determinant. We shall 
find similar singularities in their eigenfunctions. 

This work on singular Toeplitz matrices was in turn extended to give the spectrum 
of eigenvalues [H [TOj [11]. Lee, Dai, and Bettleheim [12] did the first calculation of the 
actual first order correction to the eigenvalue which they found to be of order (lnN)/N 
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when l/N was neither very close to zero or one. Their work was limited to the case 
a = 0. We shall extend their eigenvalue calculation to use in our determination of the 
eigenfunctions. 

1.3. Outline of paper 

Our study is focused on the asymptotic behavior of the eigenvalues and eigenvectors of 
the N by N Toeplitz matrices. The next chapter is devoted to finding the qualitative 
properties of the large-iV eigenvectors, calculated numerically. We note that the 
eigenstates can be classified by a momentum variable, p l , which defines the exponential 
behavior of the eigenvector. We also note that the scaling properties of the eigenvectors 
at the ends of this interval differ from the ones in the middle. We focus upon the latter. 

Chapter 3 is devoted to the justification of a kind of quasi-particle |T3], EH] theory 
for these eigenvalues and eigenvectors. Specifically we develop arguments and numerical 
data for the two equations given in the abstract. Both of these concern the "momentum" , 
p 1 . One of these defines an expression for the eigenvalue in terms of the momentum and 
N. The other says that the momenta are roughly equally spaced along the line (0, 2ir). 
Both of these are derived heuristically and checked by numerics. 

Chapter 4 makes use of the Wiener-Hopf method to construct an analytic theory 
of the eigenvectors. When a and (3 are in the right range, including for example, 
< a < —f3 < 1, we can find a closed form approximate expression for the eigenvectors 
appropriate for small j/N. In Chapter 5, this integral expression is evaluated. To carry 
out the evaluation, we must give, as input, the eigenvalue. In this chapter, we show 
that, for large N, these operator eigenvectors provide an excellent approximation to the 
corresponding matrix eigenvector for values of j/N < 0.7. Many of the scaling properties 
of the eigenvectors and eigenvalues are derived from this analysis of the operators. 

Chapter 6 briefly discusses these results and gives suggestions for future work. 

2. The eigenvectors 

2.1. Defining equations 

The Toeplitz matrices, of course, have N eigenvalues corresponding to N right 
eigenvectors, ip l - N and an equal number of left eigenvectors ■ (In discussing these 
matrices, we shall most often drop the superscript, N.) The eigenvalues and eigenvectors 



obey 





3=0 



and 



N-l 



3k 



e l ' N ftf, for all k G [0,N- 1]. 



(4b) 



3=0 
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Each of these eigenvectors are uniquely defined, up to an overall multiplicative constant, 
whenever the eigenvalues are non-degenerate. 

2.2. N = 00 

For finite- N, we do not know closed form expressions for the eigenvalues and eigenvectors 
of the Toeplitz matrix. However, the eigenvalues and eigenvectors of the Toeplitz 
matrices are very easily found when Tkj is a doubly infinite matrix with indices covering 
(— 00, 00). Because of the translational invariance in the latter situation the eigenvectors 
are of the form 

ipj = exp(ipj) and ipj = exp(-ipj). (5) 

for all possible real values of p, while the corresponding eigenvalues are determined by 
the symbol and are 

e = a(exp(— ip)). (6) 

2.3. Parity symmetry 

Note that for all finite Toeplitz matrices the right eigenvector can be calculated from 
the left eigenvector by: 

$ = <0k-i-i, foraUje[0,JV-l] (7) 

Thus, for the finite Toeplitz matrix, ijA serves in an analogous role to that played by an 
adjoint eigenvector for a Hermitian matrix. In particular, if ip 1 and ip m have different 
eigenvalues, they are orthogonal. We find that when properly normalized, the vectors 
have orthonormality relations 

N-l 

^ fttf? = 5(l,m) (8a) 

j=0 

Before normalization, the vectors are arbitrary up to multiplication by a complex 
constant. To get the normalization right one must multiply by a constant, leaving 
only an ambiguity under multiplication by ±1. The completeness condition is 

N-l 

= k) (8b) 

1=0 

We have verified Eq. flHal) and Eq. (l8bl) for the range of parameters considered in this 
paper. 

The symbol's parity transform is produced by z — > 1/z. For the Fisher-Hartwig 
symbol of Eq.Q, this change is reflected by (3 — > —f3. Thus in this case, (3 and —(3 
are equivalent for finite iV since they are obtained from one another by a symmetry 
operation. Thus they have the same eigenvalue spectrum. For infinite iV they are not 
equivalent, so that in j e [0, 00], the matrices with parameters ft and —ft have different 
eigenvalues. 
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Figure 1. Eigenvalue Distribution. The image of the symbol together with the 
eigenvalues in the complex plane for a — 1/3, [3 — —1/2; the orange curve being 
the image of the symbol; the green triangles, blue stars and purple dots being the 
eigenvalues for N — 40, 100 and 200 respectively. 



2.4- Numerical calculation of eigenvalues and eigenf unctions 

Our analysis begins by the calculation of the eigenvalues and eigenfunctions of the 
Fisher-Hartwig finite- N Toeplitz matrices. The calculations are done in Mathematica. 
When Re (a) > — |, the Toeplitz matrix elements are given by 

T k — T(j — k) = (— lV~ fe r(2a + l) f 9 ) 

Jk U ] 1 ' T(a + l + f3-j + k)T(a + l-[3 + j-k) W 

These matrix elements have the asymptotic behavior 

T(j — k) ~ — 7T- — — — , - ^ for large values of \ j — k\ 

v ' tt(j - k + (3) 2a+1 1 1 

Here the ± signs refer respectively to the cases in which j > k and j < k. 

In our numerical calculations, we use the exact form given by Eq.([9]) and then 
calculate eigenvalues and eigenfunctions using the routines supplied in the commercial 
program Mathematica. 

To set the analysis of the eigenfunctions in a proper context we first show in Figure 
[U a depiction of the eigenvalues for different N and fixed a and (3. The solid curve is 
theoretical. It shows the image of the symbol, a(z), as z traverses the unit circle. The 
dots are the eigenvalues calculated using Mathematica. According to Widom's theory 
[TUt [TT] , for many kinds of symbols, as iV goes to infinity, the spectrum of eigenvalues 
should approach that image. Specifically, the eigenvalues are approximately the ones 
given by the N = oo case in which fourier transformation gives = a(exp(— ip 1 )), 
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Figure 2. Magnitude of eigenvector. Here we plot \ipj\ against j/(N — 1) for N = 100 
and the eigenvalue with number I = [(N — l)/4]. If I ~ pN/q, with p and q being 
integers the phase goes through roughly p revolutions as j increases by q. For this 
particular value of I there are four different branches to the eigenfunction and there is 
a repetition after every four steps of j. The four branches are indicated by different 
colors and shapes in this and the next figure. The solid line is an exponential fit to 
the data. The parameter values are a = 1/3 and (3 = —1/2. 



This figure shows that, for the case plotted here, the eigenvalues follow Widom's 
prescription. In fact Eq. lfTU]) is true for all a and (3 which obey < a < |/3| < 1, 
which is the study range for this paper. 

2. 5. Different regions for eigenvector 

The next two plots, Figure [2] and Figure [3] show respectively the gross behavior of the 
magnitude and phase of these eigenvectors, one with iV = 100 and / = 24. Both plots 
show that, for large N, the ratios are constant throughout a broad central region 

of j, extending perhaps to from small values of j/(N — 1) to j/(N — 1) ~ 0.7. This 
behavior changes markedly within the first few values of j and above j/(N — 1) ~ 0.7. 
We call these two regions respectively the initial region and the final region. 



where the momenta, p l are uniformly spaced in (0,27r). In symbols, 




(10) 
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Figure 3. Phase of eigenvector. We plot the phase of the eigenvectors against 
j/(N — 1). The overall behavior is that successive points increase their phase by 
an amount close to (27r)/4. We again show the behavior as four branches, each with 
points separated by four steps of j, with the different branches being denoted by 
different colors and shapes. The branches are separated by a height difference of about 
7r/2. The phase increases by an amount slightly smaller than mir/2 as j advances by 
m units. That is why the lines slope downward. 



2.6. Central region 

The first and most obvious characteristic of the central region is that over most of the 
region the eigenvector varies exponentially 

tff « exp(w l,N j) (11) 

as indicated by the straight line behavior of exponent versus j in Figure [2] and Figure 
El This behavior can most easily be understood as a consequence of the translational 
invariance of the matrix Tj^ — T(J — k). The central region is far from the two ends 
j = and j = N — l. The eigenfunction in this region behaves as it would if we were able 
to push the two ends infinitely far away. We have picked I = [(N — l)/4], so that this 
particular eigenfunction shows four different branches which differ in phase by roughly 
tt/2. 

2. 7. Initial region 

This branch-structure can also be discerned in the initial region (See Figure HJ) In this 
plot, each branch shows a few bumps for small values of j. The bumps decrease in 
size as j increases and the branches approach one another in magnitude. (Their phases 
differ by about (27r)/4.) The behavior then settles down to the exponential decay seen 
in the central region. The behavior in the initial region does not appear to vary with N 
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Figure 4. Magnitude of eigenvector in initial region. Here we plot four branches of the 
eigenvector labeled by k = 0, 1, 2, 3. Each branch has j — fc(mod 4). For each branch, 
is plotted against j for the eigenvalue with number I = \{N — l)/4], 
N = 801, and parameter- values of a = 1/3, [3 = —1/2. The curves each show their 
different individual behavior for small j and then become almost identical to one 
another for larger j. 



for large values of N. Different values of N show ip versus j plots which have, for small 
values of j, exactly the same shape (See Figure El) 

2.8. Final region 

Figure [6] and Figure [7J show eigenfunction for relatively large values of j / {N — 1). We 
note from the former figure that the four branches have different shapes in this region. 
However, Figure [7] shows that the shapes of the different branches settle down to an iV- 
independent form as iV goes to infinity. However, the scaling in this region appears to 
be somewhat complex. We leave a further analysis of this region to a later publication. 

3. Quasi-particle results 

In the 1950s, Landau [131 [33] worked out a theory of a low-temperature Fermi liquid, 
e.g. He 3 , which argued that the Fermi liquid was just like a non- interacting Fermi gas 
except that it had a different relation between energy and momentum. Here, we would 
like to make the same kind of argument. 

3.1. "Energy" -momentum relation 

We follow Landau in saying that the relevant "momentum" , p, appears in the phase of 
the eigenfunction as exp(ipj). In the idealized case, N = oo, corresponding to Landau's 
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Figure 5. Magnitude of eigenvector in initial region for one branch. Here we plot 
ln|^j/Vj+l| 5 against j = 4fc for N = 101,201,401,801 and 1601 for the eigenvalue 
with number I = [(N — l)/4] and parameter-values of a — 1/3, (3 = —1/2. One can see 
that the shape of the curves does not change appreciably as N varies. So the initial 
region contains eigenvectors with TV-independent shapes. 




0.5 0.6 0.7 0.8 0.9 1 



j/(N-l) 

Figure 6. Magnitude of eigenvector in ending region. Here we plot ln\if}j/ipj + i\, 
against j for N — 1601 and the eigenvalue with number I — [{N— 1)/4] and parameter- 
values of a = 1/3, f3 = —1/2. The four-branch structure is clearly indicated. 
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Figure 7. Magnitude of eigenfunction in ending region for one branch. Here we plot 
In l^/Vj+lli against j = 4fc for TV = 101,201,401,801,1601 and the eigenvalue with 
number I = [(N — l)/4] and parameter-values of a — 1/3, (3 = —1/2. One can see that 
the shape of the curves does not change much as N varies. 



non-interacting case, the phase and eigenvalue obey e = a(exp(— ip)). We carry this 
result over directly to finite values of N as in Equation II of the Abstract. We write 
that equation in terms of a _1 (-), the function inverse to a(-), to define 

p{ = i\na- l (e l ) (12a) 

The momentum variable defined from the wave function will be written as p^. One 
possible definition is 

plp = -i[lnri-ln^]/(J-j) (12b) 

Here we must pick the branch of the logarithm in Eq. ffl2|) so as to make p^ approximately 
equal to 2ttI/N. We need to figure out whether, the two definitions give results that are 
substantially equal in the large 7Y limit, for example, whether it is true that 

p\=p l € + o{l/N) (13) 

Eq. (1131) and the equations in the abstract are all checked using Tables 1 through 
4. All four tables analyze the data for a = 1/3 and (3 = —1/2. The first two tables use 
data for the imaginary part of p, while the latter two refer to the real part of p. In each 
case, the p^ in the table is obtained by using Eq. (112bl) with J = [0.57Y] and j = [0.27Y]. 
The first two lines in each table describe the p-values thus obtained. They are seen to 
be increasingly close for larger values of N. The next line gives the difference between 
these values multiplied by powers of N. No one of these lines increases substantially 
as N increases. This result confirms Eq. (1131) . which states that the two methods of 
determining p are consistent with one another in the limit of large N. In fact these 
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N 


2000 


1000 


400 


200 


100 


40 


\r <p J 


0.00765381 


0.0141435 


0.0314841 


0.0570298 


0.101950 


0.213238 


9(p«) 


0.00765381 


0.0141435 


0.0314842 


0.0570311 


0.101961 


0.213463 


difference- 10 6 


-0.00111 


-0.00881 


-0.146 


-1.29 


-11.7 


-225 


difference- N 3 


-8.85 


-8.81 


-9.36 


-10.3 


-11.7 


-14.4 




15.308 


14.146 


12.594 


11.406 


10.195 


8.530 


SUM ■ N — Fit 





-0.00886 


-0.0316 


-0.0640 


-0.120 


-0.258 



Table 1. Two estimates of the imaginary part of p are shown to be substantially 
equal. They are then successfully fit to the theory. The value of I for this table is 
I = [(N — l)/2]. The third line shows the difference of the two estimates and the 
fourth line shows that difference multiplied by N 3 . The rough constancy of this line 
with increasing N supports the correctness of Eq. (fT3)) . The next to last line shows 
the result of multiplying the imaginary part of p^j, in the top row by N. According to 
the theory this quantity should be (2a + 1) In N + 9(C). We get a good fit to these 
data by using 5(C) = 2.639. The last row in the table is the result of subtracting this 
fit from the result in the previous row. It represents the error in our analysis. Since 
the relative error is so small, we argue that Equation I of the abstract is corroborated. 



N 


2000 


1000 


400 


200 


100 


40 




0.00761559 


0.0140678 


0.0313017 


0.0566988 


0.101397 


0.212772 




0.00761559 


0.0140678 


0.0313018 


0.0566981 


0.10141 


0.212644 


difference- 10 6 


0.000137 


0.00964 


-0.0699 


0.680 


-13.4 


128 


difference -N 3 


1.10 


9.64 


-4.48 


5.44 


-13.4 


8.17 


*M ■ n 


15.231 


14.068 


12.521 


11.340 


10.140 


8.511 


• JV - Fit 





-0.00811 


-0.0281 


-0.0538 


-0.0986 


-0.200 



Table 2. The same as Table [T] except that the values of I are given by I = 
[{V5-l)(N-l)/2}. 



results tend to indicate that the difference between these two estimates of p is of order 
1/iV 3 for the imaginary part and 1/N 2 for the real part. Thus we have verified Equation 
II of the abstract. 

Next we compare the p- values in the table with Equation I of the abstract. The 
next to last row of the table gives the value of N times the first row. According to 
Equation I, that row should be fitted by i(2a + l)lniV + C, where C is a complex 
constant. The final row gives the errors for each N after the corresponding theoretical 
result is subtracted away. Each table thus involves one adjustable constant. The small 
values of the remainder in the last row indicate that Equation I of the abstract appears 
to be quite justified. The fit for Tabled] is shown in graphical form in Figure El 

These results of this comparison are shown in Table [1] for p roughly equals to n/2. 
Table [2] shows the same analysis applied to an "irrational" value of p chosen by taking 
I = l(VE- 1)(JV- l)/2]. Tables [3] and a give a similar comparison for the real part of 
p. In all cases Eq.ffTBl appears to be justified. 
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Figure 8. Imaginary part of p. This plot shows the data in the next to last row of 
Table Q] plotted along with a least squares fit. The theoretical slope is 1.667 and the 
fitted one is 1.728. 



N 


2000 


1000 


400 


200 


100 


40 


ft(g^) • 10 6 


3.80031 


14.4177 


84.7912 


326.402 


1265.36 


7661.04 


R(q e ) ■ 10 6 


1.72814 


6.91045 


43.1505 


172.352 


687.673 


4274.71 


difference TO 6 


2.07 


7.51 


41.6 


154 


578 


3386 


difference -N 2 


8.29 


7.51 


6.66 


6.16 


5.78 


5.42 


ft(^) • N 


0.00760 


0.0144 


0.0339 


0.0653 


0.127 


0.306 


3fc(g^) • iV — Fit 





0.00682 


0.0263 


0.0577 


0.119 


0.299 



Table 3. The same as Table [TJ except that we examine here the real part of 
q = p — 2irl/(N — 1). The fit on the last line sets N times the real part equal to 
a constant, 0.00760. 



N 


2000 


1000 


400 


200 


100 


40 


KM ■ io 6 


-808.955 


-1610.85 


-4038.63 


-7861.35 


-14949 


-26243.8 


&{q e ) ■ 10 6 


-810.805 


-1617.64 


-3997.34 


-7726.56 


-14391.3 


-29846.7 


difference- 10 6 


1.85 


6.80 


-41.3 


-135 


-558 


3603 


difference -N 2 


7.40 


6.80 


-6.61 


-5.39 


-5.58 


5.76 


&(^) • AT 


-1.618 


-1.611 


-1.615 


-1.572 


-1.495 


-1.050 


3?(^) • N - Fit 





0.00707 


0.00246 


0.0456 


0.123 


0.568 



Table 4. The same as Table [3J except that we examine here the real part of 
q = p — 2ttI/(N — 1). The fit on the last line compares N times the real part equal to 
a constant, —1.618. As one can see, this fit is quite accurate especially for the larger 
values of N. 



CONTENTS 



14 



3.2. Spacing of eigenvalues 

The other part of quasi-particle theory that we wish to verify is that for successive l- 
values the p l should be spaced by 2n/N. This will help ensure that the phase difference 
between ip l Q and V'jv-i can remain independent of I. To check this we notice that the fits 
for the real part of p in Tables [3] and H] indicate that these two eigenvalues have their 
spacing correctly described by the two equations in the abstract. A brief examination 
shows that these spacings seem to continue to hold for all values of I except the ones 
which have l/(N — 1) <C 1 or 1 — 1/ {N — 1) <C 1. Near these end-points different spacings 
are observed. In this paper, we do not pursue this point any further. 



4. Wiener-Hopf method 

4-1. The singly infinite eigenvalue problem 

To get a good approximation for the Toeplitz matrix's eigenvector, one might start by 
looking at the limiting case in which N is infinity. The case in which the integer index, 
j, is extended from — oo to oo is simply solved by Fourier transformation as in Eq.(jSJ) 
and Eq.([6]). But the singly infinite case, described as the Toeplitz operator problem, in 
which j extends from, say, to oo is much harder. The approach from singly infinite- iV 
to the case of N large but finite is tricky and non-uniform. (We give an example of 
the non- uniformity below.) Nonetheless the singly infinite matrix provides a good entry 
into the finite- iV system. 

We want to solve the Toeplitz operator eigenvalue-problem 

oo 

'^^Tkjipj = eipt, for all k > (14a) 

j=0 

and the transpose problem 

oo 

^ipjTjk = e^k, for all k > (14b) 

j=0 

where T is a Toeplitz operator, i.e. Tjk = T(j — k). 

There is a well-developed theory, called the Wiener-Hopf method, for solving 
equations involving this kind of operator [15]. See for example [16] and Chapter IX 
of McCoy and Wu [5j. See also [TFJ [IE] for applications close to the problem studied 
here. Following McCoy and Wu, we define a class of vectors with indices which run from 
— oo to oo. The components of these special vectors are required to vanish for negative 
values of their index. The vectors of this type are indicated by the superscripts "+". 
For example the "+" vector corresponding to an eigenfunction of the Toeplitz operator 
is written as 

at = I *>• 3 a 

' I 0, j<0 
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then we can write the eigenvalue equation as 

oo 

yj Kj k ip^ = for all j > where Kj k = Tj k — e5j k 

k=—oo 

We would like to define an equation for the eigenfunction which covers the entire range, 
(— oo, oo), for the variable j. To do that we define another kind of vector, one indicated 
by a superscript "— " . One example of such a vector is 

ipj = E Kjk^k' f° r & U integer j (15) 

k=— oo 

In general, the vectors in this category vanish for j > 0. However this particular vector 
also vanishes for index equal to zero. This extra condition will play an important role 
in our solution. 

The result of this subsection is an eigenvalue equation which reads 

00 

Kj k ip£ = ip~ for all integer j (16) 

k=— oo 

This equation has the same content as Eq. (jl4ap . 
4-2. Fourier transforms 

Eq.t flBl) may be solved by Fourier transformation. Define 

00 

<z) = J2 r 0V (17) 



j=-oo 
00 



j=~oo j=0 
oo —1 

$-{z) = tj* = E ^ n ( 19 ) 

jr' = — OO j = — OC 

Multiply Eq. fll6p by and sum over all j, to give 



00 00 



E E - ^-"^ = E ^ 

j=—oo k=—oo j=—oo 

which then gives an eigenvalue equation expressed in Fourier language: 

K(z)ip + (z) = (a(z) - e) ip + (z) = tp-(z) (20) 

At this point, the eigenvalue equation has been reduced to an algebraic equation, one 
that may apparently be solved by algebraic manipulations. However, Eq. fl20l) has two 
unknowns, ip + {z) and ip~(z) so that additional knowledge is needed to obtain a solution 
for ip + (z). That knowledge is based upon the analytic properties of the various functions. 
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4-3. Factorization 

To find ip + (z), we must factorize K(z) = a(z) — e into a product of functions of the form 

where v is an integer determined by the form of K(-). It is called the winding number 
since it represents the total number of times the function wraps around zero when z is 
rotated once around zero in the positive sense. The winding number obeys 

i r , i 



2-ni 



J din K(e ip ) = ^- [In K(e 2m ) - In K (e 0i )] (22) 



o 



This factorization is important since it will enable us to immediately construct a non- 
trivial solution to the Toeplitz operator equation in the one case, v = — 1, in which this 
solution exists [16J. (See McCoy and Wu [5] page 208-215. The condition, v — — 1, 
demands that the j = component of ip~ must be zero, as is required for our solution.) 

A factorization like this can be derived simply when the function, K(z), has no 
singularities or zeros on the unit circle and has none at zero or infinity. However, all 
these conditions fail in the Fisher-Hartwig case. But, McCoy and Wu [5j argue that the 
factorization also works with the definition, Eq. (l22|) . of the winding number in the case 
in which K(z) is only continuous, not analytic, on the unit circle and has appropriate 
analyticity properties in an analytic continuation away from that circle. Their analysis 
then covers the Fisher-Hartwig case in the range interesting to us: < a < \/3\ < 1. If 
one applies this integral formula of Eq. (1221) to find the winding number for K one sees 
that 

v = -1 for - 1 < < 0, (23a) 

so this is the case in which the Toeplitz operator gives a non-trivial solution for the 
eigenvector. Conversely 

v = 1 for < < 1 (23b) 

so there are only trivial eigenvector solutions in this situation. 

McCoy and Wu also give us formulas for finding the functions K ± {z). One starts 
from the contour integration 



2iri f Sl z' — z 

Here the upper (and lower) signs respectively refer to the cases in which z is inside (and 
outside) the unit circle. Then 

K±{z) = exp[G ± (z)] (25) 

With these definitions, the eigenvalue equation has the form 

ip + (z)K + (z) = z- u ip~(z)K-(z) (26) 
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The left hand side contains only non-negative powers of z. To say when we can have 
eigenvalues, we must look at the right hand side and see whether it can match the 
left. The three factors on the right hand side contains powers which are at a maximum 
—v — 1 + 0, obtained from reading the three factors from left to right. 



44. v = -\ 

The sole possibility of obtaining a solution to Eq. (1261) comes from picking v = — 1. 
Then the right hand side can match the left by simply having each side be a constant 
independent of z. Thus one finds 

ip + (z)K + (z) = z%l)-{z)K-{z) = C (27) 

where C is a non-zero constant which will set the phase and magnitude of the wave 
function. 

In this case, we have as our solution 

^ +{z) = WGj for N < 1 (28) 

We have now obtained a solution in terms of quadratures. 

Note that in general we expect that the z/-value will remain constant as the 
eigenvalue varies over some open set. Thus one has a region in the complex plane 
in which the Toeplitz operator has eigenvalues. This result is in contrast to the outcome 
for the Toeplitz matrix which usually has N isolated eigenvalues. 



Example: A simple pole We choose a very simple example to illustrate what we have 
said. Take a to be of the Fisher-Hartwig type with a = and (3 = —1: a(z) = z~ l (See 
Eq.pi), where for convenience we have dropped the overall minus sign). The matrix, 
T jk , is 












1 







(29) 



The equation <ft = Tip is equivalent to 
4>j = if) j+1 for j > 0. 

Thus the operator T = z~ l when it is acting to the right transfers information from large 
j to smaller j. It is reasonable to assume that all operators with v — — 1 similarly transfer 
information towards smaller j, and thereby they can effectively produce a boundary 
condition of a vector being zero as j goes to infinity. 

Take the eigenvalue to be e, with magnitude less than unity, so that it sits within 
the unit circle. Then 



K(z) 



K + /(zK~ 
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The factorization then gives 



K + (z) = l-ez and K~(z) = 1. 

As required K + has its zero outside the unit circle. One therefore finds that the 
eigenfunction has a fourier transform 

C 



C (1 + ez + {ezf + ■■■) 



(30) 



1 - ez 

which in turn implies an exponentially decaying eigenfunction 

^ = C(ey for j = 0,1,2,-.. 

Thus the entire unit circle is filled with eigenvalues while ip + (z) has a single pole, and 
the eigenf unctions ipj are geometrical series. 

This situation is of the Fisher-Hartwig type with a = and (3 = — 1. 

However, the finite Fisher-Hartwig matrix with this symbol has a very different 
behavior. This matrix, has N unit entries just above the diagonal, with all the 
rest of the entries being zero. (Here the notation describes an N order Toeplitz 
matrix and not the Nth power of anything.) Thus it looks just like the infinite matrix 
of Eq. (1291) except that it is truncated after N rows and N columns, i.e. 

10 

10 

1 ■•• (31) 





T (N) 



NxN 



The eigenvalue equations read 
V>jv-i = 

and 

V>jv-i-i = (l/e)Vj\r-j for 3 = l,2...iV - 1 
These equations are solved successively. They do not make sense for e 



(32a) 



(32b) 
0. For all 

other values of e, the eigenvector has all of its elements being zero. As a result, all 
eigenfunctions are trivial. 

Thus, the eigenfunction behavior of this particular order- iV Fisher-Hartwig matrix 
matches that of the corresponding operator only in the special case in which we set 
C = in the operator eigenvector. 



5. Calculation of eigenfunctions 

The work of previous authors on eigenvalues has given us useful but limited information. 
We have an asymptotic expansion which gives the line on which the eigenvalues fall and 
their probability density on that line. Their exact placement on that line has never been 
explicitly calculated. Nor have the eigenfunctions associated with these eigenvalues. 
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-0.5 0.0 0.5 1.0 1.5 



Figure 9. Maps of the symbol on the iu-plane. The black curve is the image of the 
symbol w = a a ,f){z) for a = 1/3 and j3 = —1/2 with z going around the unit circle 
in the counterclockwise direction. The symbol has a singularity at z = 1 (mapped to 
w = 0) as we can see from the other curves in this figure. The curves counting from 
the outmost are the images of this symbol for circles \z\ = 0.75, \z\ = 0.8, \z\ — 0.85, 
|*| = 0.9, \z\ = 0.95, \z\ = 1, \z\ = 1.1, |*| = 1.2, \z\ = 1.4, \z\ = 2 and \z\ = 4 
respectively (red curves for \z\ < 1, blue curves for \z\ > 1 and the black curve for the 
unit circle). The green lines consist of the endpoints of the curves. 



Here, we set out to fill in some of these omissions, using both asymptotic calculations 
and numerical work. We shall especially make use of the Wiener-Hopf solutions outlined 
in the previous chapter to discuss the properties of the behavior arising from the Fisher- 
Hartwig generating function. 

For < a < —j3 < 1, all the eigenvalues lie within the closed curve formed as 
the image of the unit circle mapped by the symbol. Since the image winds around in 
the clockwise direction, the winding number of (a a ,p — e) is — 1 for any eigenvalue, e. 
Therefore, we can use the quadrature of the previous chapter, 

^+( z ) = Ce~ G+{z) (33) 



G+(z) = ±,<L d /^^fW ~ £)] (34) 



where 

G + (z) = o 

2m j g 4, — & 

This result is then a solution by quadratures of the Toeplitz operator problem involving 
the Fisher-Hartwig symbol. We may assume that this solution is accurate for small 
positive values of x = j/(N — 1) and becomes less accurate as x approaches one from 
below. 
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-3-2-10 1 2 3 



Figure 10. The singularities in the contour integral of in the z'-plane for the 
integration in Eq. (|34p . The original contour, on the unit circle, is shown in blue. This 
path can then be deformed into the one shown in red, encircling singularities outside 
the unit circle. The black lines are encircled branch lines, including one starting at 
z c , where a(z') = e, and another along the positive real axis. The path on the big red 
circle can be taken to infinity where it does not contribute. The dot is the singularity 
at z' = z c . There is an additional pole singularity at z' = z, which does not contribute 
because z is inside the unit circle. 



5.1. Evaluation of the integrals 

Eq.( |33l) and Eq.( l34l) together give the result for the Fisher-Hartwig symbol for the 
case < a < —f3 < 1. Figure [JJ shows the image of the unit circle formed by the 
symbol together with the eigenvalues of the corresponding Toeplitz matrices for the 
case a = 1/3, (3 = —1/2. We see that the eigenvalues are all distributed inside the 
image of the symbol. They approach that image as the size of the Toeplitz matrix goes 
to infinity. Further discussion of the asymptotic behavior of the eigenvalues has been 
given in [12] . 

To perform the contour integration in Eq.( j3~4l) . we need to rewrite our symbol in 
terms of an analytic function, which is 

a a ,p( z ) = ( z ~ l) 2Q /- a e - i(/3+a) " (35) 
We shall choose this function, which has branch points at zero, one, and infinity, to 
have a branch cut along the positive real axis (See Figure [TO]) . The logarithm in Eq. (134?) 
has an additional singularity at the point, outside the unit circle, where a(z') equals the 
eigenvalue. We call this point z c . It obeys 

a(z c ) = e (36) 

We can then have another branch line extending from z c to oo. As N goes to infinity, e 
approaches the image of the unit circle and z c approaches the unit circle. These points 
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will play important roles in what follows. Since a(z') goes exponentially to zero at 
infinity the contour in Eq. fl34|) can be deformed as shown in Figure [TUl The big circular 
path can then be taken to oo. The result is that our integral becomes, aside from an 
additive constant, 

G + {z) = L{z) + S{z) (37a) 

Here L(z) stands for a logarithmic term, which is 

L(z) = ln(z - z c ) (37b) 
while S(z) stands for a term which has a weak singularity at z — 1, 



oo 

S(z) = ± /" — In 



2m J t - z [ (t - l)**tP-<* e i(P+a)* 



(37c) 



Here, z c is the position of the zero in K(z) = a(z) — e in the z-plane. According to 
Widom's ideas, this zero approaches the unit circle as N goes to infinity, but remains 
outside of that circle. This approach makes for a near-singularity in ip + (z) on the circle. 
For z inside the circle, this near-singularity looks like a singularity just outside the point 
z = exp(— 2ml /(N — 1)). In contrast the real singularity in ip + (z) is produced by the 
branch line which passes through z — 1. The singularity and the near-singularity are 
well separated except for situations in which I is relatively close to I = or / = iV — 1. 
We have not explored these exceptional limiting cases. 



5.2. Accuracy of Wiener- H op f solution 

The solution we have generated is exact in the limiting case in which j is fixed and N 
goes to infinity. Thus, we might expect that the solution is accurate for large N and 
fixed j/(N — 1). To see how accurate it is, we plot in Figure [11] the absolute value of the 
difference between the Wiener- Hopf quadrature and the exact solution for N = 1000. 
The absolute value of the error never gets larger than a few times 10~ 7 . It decreases 
slowly as j/(N — 1) increases. The relative error, plotted in Figure [121 grows with 
j/(N — 1), starting from order 10 -7 and increasing to order one at j/(N — 1) ~ 1. 
Therefore, we may say that the Wiener-Hopf solution reproduces the main features of 
the exact eigenfunction except very close to the maximum values of j/(N — 1). 

5. 3. Examination of solution 

Even a superficial examination of Eq. (1371) shows that its behavior fits the known 
properties of the eigenfunction. One term in the wave function is 

1 1 00 

Mz) = exp (-L(z)) = = — E 0*/^ ( 38 ) 

Z C {1 Z/Z c ) z c _. =Q 

This part of ip{z) thus produces a ipj which is 

(^ = -(l/^ +1 (39) 
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Figure 11. Absolute comparison of Wiener-Hopf solution with exact eigenfunction. 
The magnitude of the difference between the Wiener-Hopf solution and the exact 
eigenfunction is plotted against j/(N - 1) for N = 1000, I — [(N — l)/4], a — 1/3 and 
P = -1/2. 
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Figure 12. Relative comparison of Wiener-Hopf solution with exact eigenfunction. 
The relative magnitude, \ip/ip WH — l| of the difference between the Wiener-Hopf 
solution and the exact eigenfunction is plotted against j/(N — 1) for N = 1000, 
I = [(N — l)/4], a = 1/3 and (3 = -1/2. 
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Figure 13. Expansion of rpg. The magnitude of the terms in the expansion of ips{z) 
in powers of z. These terms reflect corrections to the eigenfunction beyond the leading 
exponential behavior. The notable features are a long tail proportional to (J + l) _2a_1 
and a few decreasing order-one bumps for small j. 



This structure is exactly the decaying exponential so evident in the central region of j. 

The other parts of ip(z) come from ips{z) — exp(—S(z)). This exponential is 
independent of N except for the very weak iV-dependence produced by a very small 
additive term of the form j(hiN)/N in the eigenvalue, e. This effect of the N- 
dependence in this term is quite negligible. 

Figure [T3l plots the magnitude of the terms in the power series expansion of il>s( z ) 
in z. Using the same notation as used for 4>l(z), these terms are called (ips)j- The first 
few terms are of order unity and produce the small- j bumps evident in Figure HI The 
bumps slowly die away for higher j. 

For higher j the behavior of (ips)j is dominated by a power law behavior, which 
shows up as a straight line in Figure [13j Even though these terms fall off, this 
contribution remains significant. This behavior is shown in Eq.( )40l) . It results from 
a singularity as z goes to unity evidenced in Eq. fl37cl) , which shows a behavior for small 
z — 1 of the form 

lnfeM _^(°+fl [ L ds ^ 



IT £ Jo 1 + S — Z 

L is a large number which serves as a cutoff at large values of s. This expression can be 
expanded in a power series in z, with coefficients of high order terms in the form 

smTr(a + (3) f L s 2a 



sinvr (a + (3) f 
J 71 e Jo 



For large j the integral contributes for small values of s, of order 1/j so that we can 
replace (1 + s)^ 1 by e s ^ +l > . We thus get the leading singularity in h\ips{ z ) to be given 
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by 

n , / m T(2a + 1) sin7r(a + (3) , inS 

y*M*))i - - 7re( ; + 1)^! ( 4 °) 

Since this singularity in ]ni[)g(z) indicates a zero added onto a leading term which is a 
constant, the very same form of singularity appears in ipg(z), which therefore has a high 
order expansion 

T(2a + l)sin7r(a + /5) 
IM'h ne{j + 1)2a+1 (41) 

for large values of j. 
5-4- Mixed terms 

5.4-1- Sums So far we have argued that the solution for ipj should contain terms of the 
exponential form shown in Eq. (1391) as well as the algebraic form in Eq.f HTT) . In addition 
there is a set of bumps which appear only for small j. For larger j, then, we might 
expect an eigenfunction which looks like 

i>j ~ A{z c )~ i ~ 1 + B(j + l)-( 2a+1 ) (42) 

where A and B are both constants. Since both terms are generated near j = 0, we 
expect both A and B to be of order unity. Note that both terms fall off for higher j. If, 
as we expect, 3?(lnl/^ c ) is small, then the exponential decay will be slow for small j, 
while the decay of the algebraic term will be rapid. However, for sufficiently large j the 
two terms return to being of the same order. We assume this equality is achieved when 
j = O(N) so that we might have a boundary condition in which the two terms effectively 
interfere to produce a very small result near and beyond the border at j — N — 1: 

(z c y N ~ aH2«+d 

Since z c = exp(— ip), we can calculate the imaginary part of p as 

3 (p) = ft (ln* c ) = (2a + 1)^ + 0(1/N) (43) 

The 0(1/N) term arises from the ratio of A/B in Eq.( ]42l . This change in p enables 
one to calculate the shift in the eigenvalue in the form 

5e dine 

— = -k— $ iP) 
e op 

so that the change in the eigenvalue is given by 

-=^P + ia cot |J (2a + 1)— (44) 

This is, to our knowledge, the first calculation of the eigenvalue shift for a not equal to 
zero. 
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5.4-2. Products The actual form of the Wiener-Hopf eigenfunction is not a sum of the 
form shown in Eq. fj42l but instead a product in which 

1>(z) ~{z- zj-^siz) (45) 

which can be Fourier transformed to give 

„ T(2a + l)sinn(a + ( 3) > 1 

We introduce an integral representation for the last factor to give 

sin7r(a + (3) 
re e z, 



JO ,_ n 



^ -(k+i)ij. v -j+k 

fc=0 

The geometric series can be summed to find 
sin 7i(a + (3) 



dfi e -> 2a ^ ! 



1 - 2 c e-^ 

This elegant result enables us to write ipj as a sum of two terms: 

^ ~ Az-J- 1 + B(j + l)-( 2a+1 ) (47a) 

where the term multiplied by A dominates in the central region and both terms are of 
the same order in the ending region. Here the value of the coefficients are 
sin7r(a + /3) f°° J __ _ 2a 1 



A 1 / / dfi e~^ 2a - r — - (47b) 

and 



7r e ./n 1 — ^ c e M 



sinvr a + f , „ 2n 1 
£ — / ^ e"> 2 a - j-— r 47c 

To construct this form for B we have made the change of variables \x — > + 1). 
Note that this expression gives an eigenfunction of precisely the same form as we have 
analyzed in Eq. (1421) in the previous subsection. Therefore the previous analysis will 
hold, if we assume that the coefficients A and B are constants. The only problem is 
that the integrals for A and B have singularities if z c is close to unity. For large iV this 
value of z c will only occur if l/(N — 1) or (1 — l/(N — 1)) is very small. If we are not 
in those special regions of I, then A and B can be considered to be constants of order 
unity so that our previous estimate of the eigenvalue given by Eq.fflJ will hold. 



6. More to investigate 

This paper is only a partial investigation of the problem at hand. It focuses on obtaining 
a good solution for smaller values of j/(N — 1) and does not look in detail at values of 
j/(N — 1) very close to one. In particular, the form of large-j scaling is not considered. 
One might expect that some sort of Wiener-Hopf solution might be obtained on the 
large- j side and that values of p l and e l might be calculable through order 1/JV. 

No calculations are extended to the case a = 0, nor to the case in which \/3\ < a. 
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We have shed no light on the source of the Fisher-Hartwig a 2 —(3 2 in the determinant 
of T nor have we explained what happens when l/(N — 1) is close to zero or one. 

Nonetheless this paper does contain a large number of new results: we see excellent 
eigenfunctions for smaller j, and corrections in p l and e l to order (In N)/N. More 
importantly we have shed some light on how the logarithms arise to permit the 
interference between an exponential and an algebraic term in the eigenfunction. In 
all of this, we have made a novel use of "quasi-particle" theory a la Landau. 
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